library(texreg)
library(dplyr)
library(doBy)
library(lme4)

the_prefix <- ""

mean.na <- function(x){
  mean(x,na.rm=T)
}


### geographically merged file; generated by:
### table3-01-HIX-merge-12112020.R
load("kff-allgeocodes-12212018.Rdata")

### recode variables
gdta.co$AGEN <- as.numeric(as.character(gdta.co$AGE))

gdta.co$mean_premi27.15.14.D <- gdta.co$mean_premi27.15-gdta.co$mean_premi27.14
gdta.co$mean_premi27.16.15.D <- gdta.co$mean_premi27.16-gdta.co$mean_premi27.15
gdta.co$mean_premi27.17.16.D <- gdta.co$mean_premi27.17-gdta.co$mean_premi27.16

gdta.co.15 <- gdta.co[gdta.co$DATE %in% c("JAN15","FEB15","MAR15","APR15","MAY15","JUN15","JUL15","AUG15","SEP15","OCT15","NOV15","DEC15"),]
gdta.co.16 <- gdta.co[gdta.co$DATE %in% c("JAN16","FEB16","MAR16","APR16","MAY16","JUN16","JUL16","AUG16","SEP16","OCT16","NOV16","DEC16"),]
gdta.co.17 <- gdta.co[gdta.co$DATE %in% c("JAN17","FEB17","MAR17","APR17","MAY17","JUN17","JUL17","AUG17","SEP17","OCT17","NOV17","DEC17"),]

names.11 <- c("JAN11","FEB11","MAR11","APR11","MAY11","JUN11","JUL11","AUG11","SEP11","OCT11","NOV11","DEC11")
names.12 <- c("JAN12","FEB12","MAR12","APR12","MAY12","JUN12","JUL12","AUG12","SEP12","OCT12","NOV12","DEC12")
names.13 <- c("JAN13","FEB13","MAR13","APR13","MAY13","JUN13","JUL13","AUG13","SEP13","OCT13","NOV13","DEC13")
names.14 <- c("JAN14","FEB14","MAR14","APR14","MAY14","JUN14","JUL14","AUG14","SEP14","OCT14","NOV14","DEC14")
names.15 <- c("JAN15","FEB15","MAR15","APR15","MAY15","JUN15","JUL15","AUG15","SEP15","OCT15","NOV15","DEC15")
names.16 <- c("JAN16","FEB16","MAR16","APR16","MAY16","JUN16","JUL16","AUG16","SEP16","OCT16","NOV16","DEC16")
names.17 <- c("JAN17","FEB17","MAR17","APR17","MAY17","JUN17","JUL17","AUG17","SEP17","OCT17","NOV17","DEC17")

### descriptive statistics
vars.15 <- c("count_carrier.15","count_plan.15","count_gold.15","count_silver.15","count_bronze.15",       "count_platinum.15","count_catastr.15","count_ppo.15","count_hmo.15","count_epo.15","count_pos.15","min_premi27.15","max_premi27.15","mean_premi27.15","premi27_bench.15","premi27_sil_lowest.15","mean_premi27.15.14.D")
vars.16 <- c("count_carrier.16","count_plan.16","count_gold.16","count_silver.16","count_bronze.16",       "count_platinum.16","count_catastr.16","count_ppo.16","count_hmo.16","count_epo.16","count_pos.16","min_premi27.16","max_premi27.16","mean_premi27.16","premi27_bench.16","premi27_sil_lowest.16","mean_premi27.16.15.D")
vars.17 <- c("count_carrier.17","count_plan.17","count_gold.17","count_silver.17","count_bronze.17",       "count_platinum.17","count_catastr.17","count_ppo.17","count_hmo.17","count_epo.17","count_pos.17","min_premi27.17","max_premi27.17","mean_premi27.17","premi27_bench.17","premi27_sil_lowest.17","mean_premi27.17.16.D")


mean.vec.15 <- apply(gdta.co.15[,vars.15],2,mean.na)
mean.vec.15
mean.vec.16 <- apply(gdta.co.16[,vars.16],2,mean.na)
mean.vec.16
mean.vec.17 <- apply(gdta.co.17[,vars.17],2,mean.na)
mean.vec.17

gdta.co$YEAR <- NA
gdta.co$YEAR[gdta.co$DATE %in% names.11] <- 2011
gdta.co$YEAR[gdta.co$DATE %in% names.12] <- 2012
gdta.co$YEAR[gdta.co$DATE %in% names.13] <- 2013
gdta.co$YEAR[gdta.co$DATE %in% names.14] <- 2014
gdta.co$YEAR[gdta.co$DATE %in% names.15] <- 2015
gdta.co$YEAR[gdta.co$DATE %in% names.16] <- 2016
gdta.co$YEAR[gdta.co$DATE %in% names.17] <- 2017

mean.na <- function(x){
	mean(x,na.rm=T)
}

gdta.co$high.mean.14 <- 1*(gdta.co$mean_premi27.14 > quantile(gdta.co$mean_premi27.14,na.rm=T,.75))
gdta.co$high.mean.15 <- 1*(gdta.co$mean_premi27.15 > quantile(gdta.co$mean_premi27.15,na.rm=T,.75))
gdta.co$high.d.15 <- 1*(gdta.co$mean_premi27.15.14.D > quantile(gdta.co$mean_premi27.15.14.D,na.rm=T,.75))

gdta.co$high.mean.16 <- 1*(gdta.co$mean_premi27.16 > quantile(gdta.co$mean_premi27.16,na.rm=T,.75))
gdta.co$high.d.16 <- 1*(gdta.co$mean_premi27.16.15.D > quantile(gdta.co$mean_premi27.16.15.D,na.rm=T,.75))

gdta.co$high.mean.17 <- 1*(gdta.co$mean_premi27.17 > quantile(gdta.co$mean_premi27.17,na.rm=T,.75))
gdta.co$high.d.17 <- 1*(gdta.co$mean_premi27.17.16.D > quantile(gdta.co$mean_premi27.17.16.D,na.rm=T,.75))

gdta.co$high.year <- 0
gdta.co$high.year[gdta.co$YEAR==2017 & gdta.co$high.mean.17==1] <- 1
gdta.co$high.year[gdta.co$YEAR==2016 & gdta.co$high.mean.16==1] <- 1
gdta.co$high.year[gdta.co$YEAR==2015 & gdta.co$high.mean.15==1] <- 1
gdta.co$high.year[gdta.co$YEAR==2014 & gdta.co$high.mean.14==1] <- 1
gdta.co$high.year[gdta.co$YEAR==2013 & gdta.co$high.mean.14==1] <- 1
gdta.co$high.year[gdta.co$YEAR==2012 & gdta.co$high.mean.14==1] <- 1
gdta.co$high.year[gdta.co$YEAR==2011 & gdta.co$high.mean.14==1] <- 1

gdta.co.sub <- gdta.co[gdta.co$YEAR %in% c(2011:2017),]

### statistical models
### 2015
### model for non-market
lout15 <- lm(FAVOR ~ EDUC+as.factor(PID5)+INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+I(COPCUNEMP16-COPCUNEMP10)+COPCBCH10+as.factor(DATN)+count_plan.15+mean_premi27.15+mean_premi27.15.14.D,data=gdta.co.15[gdta.co.15$MARKET==0,])

### model for market respondents
dta.15.mkt <- gdta.co.15[gdta.co.15$MARKET==1,]
lout15a <- lm(FAVOR ~ EDUC+as.factor(PID5)+INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+I(COPCUNEMP16-COPCUNEMP10)+COPCBCH10+as.factor(DATN)+count_plan.15+mean_premi27.15+mean_premi27.15.14.D,data=gdta.co.15[gdta.co.15$MARKET==1,])

### model for uninsured
lout15u <- lm(FAVOR ~ EDUC+as.factor(PID5)+INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+I(COPCUNEMP16-COPCUNEMP10)+COPCBCH10+as.factor(DATN)+count_plan.15+mean_premi27.15+mean_premi27.15.14.D,data=gdta.co.15[gdta.co.15$COVERED==0,])

### extract data for subsets for 2015
mf15 <- model.frame(lout15)
mf15a <- model.frame(lout15a)
mf15u <- model.frame(lout15u)

### 2016
lout16 <- lm(FAVOR ~
EDUC+as.factor(PID5)+INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+I(COPCUNEMP16-COPCUNEMP10)+COPCBCH10+as.factor(DATN)+count_plan.16+mean_premi27.16+mean_premi27.16.15.D,data=gdta.co.16[gdta.co.16$MARKET==0,])

lout16a <- lm(FAVOR ~
EDUC+as.factor(PID5)+INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+I(COPCUNEMP16-COPCUNEMP10)+COPCBCH10+as.factor(DATN)+count_plan.16+mean_premi27.16+mean_premi27.16.15.D,data=gdta.co.16[gdta.co.16$MARKET==1,])

lout16u <- lm(FAVOR ~
EDUC+as.factor(PID5)+INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+I(COPCUNEMP16-COPCUNEMP10)+COPCBCH10+as.factor(DATN)+count_plan.16+mean_premi27.16+mean_premi27.16.15.D,data=gdta.co.16[gdta.co.16$COVERED==0,])

mf16 <- model.frame(lout16)
mf16a <- model.frame(lout16a)
mf16u <- model.frame(lout16u)

lout17 <- lm(FAVOR ~
EDUC+as.factor(PID5)+INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+I(COPCUNEMP16-COPCUNEMP10)+COPCBCH10+as.factor(DATN)+count_plan.17+mean_premi27.17+mean_premi27.17.16.D,data=gdta.co.17[gdta.co.17$MARKET==0,])

lout17a <- lm(FAVOR ~
EDUC+as.factor(PID5)+INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+I(COPCUNEMP16-COPCUNEMP10)+COPCBCH10+as.factor(DATN)+count_plan.17+mean_premi27.17+mean_premi27.17.16.D,data=gdta.co.17[gdta.co.17$MARKET==1,])

lout17u <- lm(FAVOR ~
EDUC+as.factor(PID5)+INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+I(COPCUNEMP16-COPCUNEMP10)+COPCBCH10+as.factor(DATN)+count_plan.17+mean_premi27.17+mean_premi27.17.16.D,data=gdta.co.17[gdta.co.17$COVERED==0,])

mf17 <- model.frame(lout17)
mf17a <- model.frame(lout17a)
mf17u <- model.frame(lout17u)

### joint
mf15$YEAR <- 2015
mf16$YEAR <- 2016
mf17$YEAR <- 2017

mf15a$YEAR <- 2015
mf16a$YEAR <- 2016
mf17a$YEAR <- 2017

mf15u$YEAR <- 2015
mf16u$YEAR <- 2016
mf17u$YEAR <- 2017

colnames(mf17) <- colnames(mf16) <- colnames(mf15) <- c("FAVOR","EDUC","PID5","INCOME","BLACK","HISP","ASIAN","AGEN","COPCNHBL10","COPCHIS10",
                                                        "COMHINC10","COPCPOOR10","COPCUNEMP10","COPCUNEMPD1016","COPCBCH10","DATN",
                                                        "count_plan","mean_premi27","d_mean_premi27","YEAR")
colnames(mf17a) <- colnames(mf16a) <- colnames(mf15a) <- c("FAVOR","EDUC","PID5","INCOME","BLACK","HISP","ASIAN","AGEN","COPCNHBL10","COPCHIS10",
                                                           "COMHINC10","COPCPOOR10","COPCUNEMP10","COPCUNEMPD1016","COPCBCH10","DATN",
                                                           "count_plan","mean_premi27","d_mean_premi27","YEAR")
colnames(mf17u) <- colnames(mf16u) <- colnames(mf15u) <- c("FAVOR","EDUC","PID5","INCOME","BLACK","HISP","ASIAN","AGEN","COPCNHBL10","COPCHIS10",
                                                           "COMHINC10","COPCPOOR10","COPCUNEMP10","COPCUNEMPD1016","COPCBCH10","DATN",
                                                           "count_plan","mean_premi27","d_mean_premi27","YEAR")

### join data by rows
mfj.u <- rbind(mf15u,mf16u,mf17u)
mfj.nmkt <- rbind(mf15,mf16,mf17)
mfj <- rbind(mf15a,mf16a,mf17a)

#### estimate models for joint data
mjf.all <- rbind(
    mfj.u %>% mutate(MARKET=NA, COVERED=0),
    mfj.nmkt %>% mutate(MARKET=0, COVERED=NA),
    mfj %>% mutate(MARKET=1, COVERED=NA)
)

lout <- lmer(
    FAVOR ~as.factor(DATN)+as.factor(PID5)+EDUC+
INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+COPCUNEMPD1016+COPCBCH10+
scale_log_count_plan+scale_mean_premi27+scale_d_mean_premi27 +(1|YEAR),
data=mjf.all %>% mutate(
                     scale_log_count_plan=as.vector(scale(log(count_plan+1))),
                     scale_mean_premi27=as.vector(scale(mean_premi27)),
                     scale_d_mean_premi27=as.vector(scale(d_mean_premi27))
                 ) %>%
filter(MARKET==1 & is.na(COVERED)),
)

lout.u <- lmer(
    FAVOR ~as.factor(DATN)+as.factor(PID5)+EDUC+
INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+COPCUNEMPD1016+COPCBCH10+
scale_log_count_plan+scale_mean_premi27+scale_d_mean_premi27 +(1|YEAR),
data=mjf.all %>% mutate(
                     scale_log_count_plan=as.vector(scale(log(count_plan+1))),
                     scale_mean_premi27=as.vector(scale(mean_premi27)),
                     scale_d_mean_premi27=as.vector(scale(d_mean_premi27))
                 ) %>%
filter(COVERED==0 & is.na(MARKET)),
)

lout.nomkt <- lmer(
    FAVOR ~as.factor(DATN)+as.factor(PID5)+EDUC+
INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+COPCUNEMPD1016+COPCBCH10+
scale_log_count_plan+scale_mean_premi27+scale_d_mean_premi27 +(1|YEAR),
data=mjf.all %>% mutate(
                     scale_log_count_plan=as.vector(scale(log(count_plan+1))),
                     scale_mean_premi27=as.vector(scale(mean_premi27)),
                     scale_d_mean_premi27=as.vector(scale(d_mean_premi27))
                 ) %>%
filter(MARKET==0 & is.na(COVERED)),
)

library(stargazer)
## stargazer(
##     lout.nomkt,lout,lout.u,
##     digits=2,
##     star.cutoffs = c(0.05, 0.01, 0.001),## ,star.cutoffs=,
##     omit=paste(c(
##         "DATN","EDUC","INC","BLACK","HISP","ASIAN","AGEN","CO","PID",
##         "Constant"
##     ), collapse="|"),
##     covariate.labels=c(
##         "\\multirow{2}{4cm}{Number of plans (logged, in sd's)}",
##         "\\multirow{2}{4cm}{Mean Premium (in sd's)",
##         "\\multirow{2}{4cm}{Mean Change in Premium (in sd's)}"
##     ),
##     omit.stat=c("ser","f","adj.rsq", "rsq"),
##     table.layout = "ts", float=F
## )

## stargazer(
##     lout.nomkt,lout,lout.u,
##     digits=2,
##     star.cutoffs = c(0.05, 0.01, 0.001),## ,star.cutoffs=,
##     ## omit=paste(c(
##     ##     "DATN","EDUC","INC","BLACK","HISP","ASIAN","AGEN","CO","PID",
##     ##     "Constant"
##     ## ), collapse="|"),
##     omit.stat=c("ser","f","adj.rsq", "rsq"),
##     table.layout = "ts", float=F
## )

### re-run models
### including Benchmark Plan data ----

gdta.co$premi27_bench.15.14.D <- gdta.co$premi27_bench.15-gdta.co$premi27_bench.14
gdta.co$premi27_bench.16.15.D <- gdta.co$premi27_bench.16-gdta.co$premi27_bench.15
gdta.co$premi27_bench.17.16.D <- gdta.co$premi27_bench.17-gdta.co$premi27_bench.16

gdta.co.15 <- gdta.co[gdta.co$DATE %in% c("JAN15","FEB15","MAR15","APR15","MAY15","JUN15","JUL15","AUG15","SEP15","OCT15","NOV15","DEC15"),]
gdta.co.16 <- gdta.co[gdta.co$DATE %in% c("JAN16","FEB16","MAR16","APR16","MAY16","JUN16","JUL16","AUG16","SEP16","OCT16","NOV16","DEC16"),]
gdta.co.17 <- gdta.co[gdta.co$DATE %in% c("JAN17","FEB17","MAR17","APR17","MAY17","JUN17","JUL17","AUG17","SEP17","OCT17","NOV17","DEC17"),]

names.11 <- c("JAN11","FEB11","MAR11","APR11","MAY11","JUN11","JUL11","AUG11","SEP11","OCT11","NOV11","DEC11")
names.12 <- c("JAN12","FEB12","MAR12","APR12","MAY12","JUN12","JUL12","AUG12","SEP12","OCT12","NOV12","DEC12")
names.13 <- c("JAN13","FEB13","MAR13","APR13","MAY13","JUN13","JUL13","AUG13","SEP13","OCT13","NOV13","DEC13")
names.14 <- c("JAN14","FEB14","MAR14","APR14","MAY14","JUN14","JUL14","AUG14","SEP14","OCT14","NOV14","DEC14")
names.15 <- c("JAN15","FEB15","MAR15","APR15","MAY15","JUN15","JUL15","AUG15","SEP15","OCT15","NOV15","DEC15")
names.16 <- c("JAN16","FEB16","MAR16","APR16","MAY16","JUN16","JUL16","AUG16","SEP16","OCT16","NOV16","DEC16")
names.17 <- c("JAN17","FEB17","MAR17","APR17","MAY17","JUN17","JUL17","AUG17","SEP17","OCT17","NOV17","DEC17")

### descriptive statistics
vars.15 <- c("count_carrier.15","count_plan.15","count_gold.15","count_silver.15","count_bronze.15",       "count_platinum.15","count_catastr.15","count_ppo.15","count_hmo.15","count_epo.15","count_pos.15","min_premi27.15","max_premi27.15","mean_premi27.15","premi27_bench.15","premi27_sil_lowest.15","mean_premi27.15.14.D")
vars.16 <- c("count_carrier.16","count_plan.16","count_gold.16","count_silver.16","count_bronze.16",       "count_platinum.16","count_catastr.16","count_ppo.16","count_hmo.16","count_epo.16","count_pos.16","min_premi27.16","max_premi27.16","mean_premi27.16","premi27_bench.16","premi27_sil_lowest.16","mean_premi27.16.15.D")
vars.17 <- c("count_carrier.17","count_plan.17","count_gold.17","count_silver.17","count_bronze.17",       "count_platinum.17","count_catastr.17","count_ppo.17","count_hmo.17","count_epo.17","count_pos.17","min_premi27.17","max_premi27.17","mean_premi27.17","premi27_bench.17","premi27_sil_lowest.17","mean_premi27.17.16.D")

mean.vec.15 <- apply(gdta.co.15[,vars.15],2,mean.na)
mean.vec.15
mean.vec.16 <- apply(gdta.co.16[,vars.16],2,mean.na)
mean.vec.16
mean.vec.17 <- apply(gdta.co.17[,vars.17],2,mean.na)
mean.vec.17

gdta.co$YEAR <- NA
gdta.co$YEAR[gdta.co$DATE %in% names.11] <- 2011
gdta.co$YEAR[gdta.co$DATE %in% names.12] <- 2012
gdta.co$YEAR[gdta.co$DATE %in% names.13] <- 2013
gdta.co$YEAR[gdta.co$DATE %in% names.14] <- 2014
gdta.co$YEAR[gdta.co$DATE %in% names.15] <- 2015
gdta.co$YEAR[gdta.co$DATE %in% names.16] <- 2016
gdta.co$YEAR[gdta.co$DATE %in% names.17] <- 2017

gdta.co$high.mean.14 <- 1*(gdta.co$mean_premi27.14 > quantile(gdta.co$mean_premi27.14,na.rm=T,.75))
gdta.co$high.mean.15 <- 1*(gdta.co$mean_premi27.15 > quantile(gdta.co$mean_premi27.15,na.rm=T,.75))
gdta.co$high.d.15 <- 1*(gdta.co$mean_premi27.15.14.D > quantile(gdta.co$mean_premi27.15.14.D,na.rm=T,.75))

gdta.co$high.mean.16 <- 1*(gdta.co$mean_premi27.16 > quantile(gdta.co$mean_premi27.16,na.rm=T,.75))
gdta.co$high.d.16 <- 1*(gdta.co$mean_premi27.16.15.D > quantile(gdta.co$mean_premi27.16.15.D,na.rm=T,.75))

gdta.co$high.mean.17 <- 1*(gdta.co$mean_premi27.17 > quantile(gdta.co$mean_premi27.17,na.rm=T,.75))
gdta.co$high.d.17 <- 1*(gdta.co$mean_premi27.17.16.D > quantile(gdta.co$mean_premi27.17.16.D,na.rm=T,.75))

gdta.co$high.year <- 0
gdta.co$high.year[gdta.co$YEAR==2017 & gdta.co$high.mean.17==1] <- 1
gdta.co$high.year[gdta.co$YEAR==2016 & gdta.co$high.mean.16==1] <- 1
gdta.co$high.year[gdta.co$YEAR==2015 & gdta.co$high.mean.15==1] <- 1
gdta.co$high.year[gdta.co$YEAR==2014 & gdta.co$high.mean.14==1] <- 1
gdta.co$high.year[gdta.co$YEAR==2013 & gdta.co$high.mean.14==1] <- 1
gdta.co$high.year[gdta.co$YEAR==2012 & gdta.co$high.mean.14==1] <- 1
gdta.co$high.year[gdta.co$YEAR==2011 & gdta.co$high.mean.14==1] <- 1

gdta.co.sub <- gdta.co[gdta.co$YEAR %in% c(2011:2017),]

### model for non-market
lout15 <- lm(FAVOR ~ EDUC+as.factor(PID5)+INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+I(COPCUNEMP16-COPCUNEMP10)+COPCBCH10+as.factor(DATN)+count_plan.15+mean_premi27.15+mean_premi27.15.14.D+premi27_bench.15.14.D+premi27_bench.15,data=gdta.co.15[gdta.co.15$MARKET==0,])

### model for market respondents
dta.15.mkt <- gdta.co.15[gdta.co.15$MARKET==1,]
lout15a <- lm(FAVOR ~ EDUC+as.factor(PID5)+INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+I(COPCUNEMP16-COPCUNEMP10)+COPCBCH10+as.factor(DATN)+count_plan.15+mean_premi27.15+mean_premi27.15.14.D+premi27_bench.15.14.D+premi27_bench.15,data=gdta.co.15[gdta.co.15$MARKET==1,])

### model for uninsurer
lout15u <- lm(FAVOR ~ EDUC+as.factor(PID5)+INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+I(COPCUNEMP16-COPCUNEMP10)+COPCBCH10+as.factor(DATN)+count_plan.15+mean_premi27.15+mean_premi27.15.14.D+premi27_bench.15.14.D+premi27_bench.15,data=gdta.co.15[gdta.co.15$COVERED==0,])
texreg(list(lout15,lout15a,lout15u),digits=4,stars=0.05)

mf15 <- model.frame(lout15)
mf15a <- model.frame(lout15a)
mf15u <- model.frame(lout15u)

gdta.co.15.mkt <- gdta.co.15[gdta.co.15$MARKET==1 & ! gdta.co.15$MARKET %in% c(NA),]

lout16 <- lm(FAVOR ~
               EDUC+as.factor(PID5)+INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+I(COPCUNEMP16-COPCUNEMP10)+COPCBCH10+as.factor(DATN)+count_plan.16+mean_premi27.16+mean_premi27.16.15.D+premi27_bench.16.15.D+premi27_bench.16,
             data=gdta.co.16[gdta.co.16$MARKET==0,])

lout16a <- lm(FAVOR ~
                EDUC+as.factor(PID5)+INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+I(COPCUNEMP16-COPCUNEMP10)+COPCBCH10+as.factor(DATN)+count_plan.16+mean_premi27.16+mean_premi27.16.15.D+premi27_bench.16.15.D+premi27_bench.16,
              data=gdta.co.16[gdta.co.16$MARKET==1,])

lout16u <- lm(FAVOR ~
                EDUC+as.factor(PID5)+INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+I(COPCUNEMP16-COPCUNEMP10)+COPCBCH10+as.factor(DATN)+count_plan.16+mean_premi27.16+mean_premi27.16.15.D+premi27_bench.16.15.D+premi27_bench.16,
              data=gdta.co.16[gdta.co.16$COVERED==0,])

mf16 <- model.frame(lout16)
mf16a <- model.frame(lout16a)
mf16u <- model.frame(lout16u)

lout17 <- lm(FAVOR ~
               EDUC+as.factor(PID5)+INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+I(COPCUNEMP16-COPCUNEMP10)+COPCBCH10+
               as.factor(DATN)+count_plan.17+mean_premi27.17+mean_premi27.17.16.D+premi27_bench.17.16.D+premi27_bench.17,data=gdta.co.17[gdta.co.17$MARKET==0,])

lout17a <- lm(FAVOR ~
                EDUC+as.factor(PID5)+INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+I(COPCUNEMP16-COPCUNEMP10)+COPCBCH10+
                as.factor(DATN)+count_plan.17+mean_premi27.17+mean_premi27.17.16.D+premi27_bench.17.16.D+premi27_bench.17,data=gdta.co.17[gdta.co.17$MARKET==1,])

lout17u <- lm(FAVOR ~
                EDUC+as.factor(PID5)+INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+I(COPCUNEMP16-COPCUNEMP10)+COPCBCH10+
                as.factor(DATN)+count_plan.17+mean_premi27.17+mean_premi27.17.16.D+premi27_bench.17.16.D++premi27_bench.17,data=gdta.co.17[gdta.co.17$COVERED==0,])

mf17 <- model.frame(lout17)
mf17a <- model.frame(lout17a)
mf17u <- model.frame(lout17u)

### joint
mf15$YEAR <- 2015
mf16$YEAR <- 2016
mf17$YEAR <- 2017

mf15a$YEAR <- 2015
mf16a$YEAR <- 2016
mf17a$YEAR <- 2017

mf15u$YEAR <- 2015
mf16u$YEAR <- 2016
mf17u$YEAR <- 2017

cn.test <-  c("FAVOR","EDUC","PID5","INCOME","BLACK","HISP","ASIAN","AGEN","COPCNHBL10","COPCHIS10",
              "COMHINC10","COPCPOOR10","COPCUNEMP10","COPCUNEMPD1016","COPCBCH10","DATN",
              "count_plan","mean_premi27","d_mean_premi27","d_bench_premi27","premi27_bench","YEAR")

cbind(cn.test,colnames(mf17u))

colnames(mf17) <- colnames(mf16) <- colnames(mf15) <- c("FAVOR","EDUC","PID5","INCOME","BLACK","HISP","ASIAN","AGEN","COPCNHBL10","COPCHIS10",
                                                        "COMHINC10","COPCPOOR10","COPCUNEMP10","COPCUNEMPD1016","COPCBCH10","DATN",
                                                        "count_plan","mean_premi27","d_mean_premi27","d_bench_premi27","premi27_bench","YEAR")
colnames(mf17a) <- colnames(mf16a) <- colnames(mf15a) <- c("FAVOR","EDUC","PID5","INCOME","BLACK","HISP","ASIAN","AGEN","COPCNHBL10","COPCHIS10",
                                                           "COMHINC10","COPCPOOR10","COPCUNEMP10","COPCUNEMPD1016","COPCBCH10","DATN",
                                                           "count_plan","mean_premi27","d_mean_premi27","d_bench_premi27","premi27_bench","YEAR")
colnames(mf17u) <- colnames(mf16u) <- colnames(mf15u) <- c("FAVOR","EDUC","PID5","INCOME","BLACK","HISP","ASIAN","AGEN","COPCNHBL10","COPCHIS10",
                                                           "COMHINC10","COPCPOOR10","COPCUNEMP10","COPCUNEMPD1016","COPCBCH10","DATN",
                                                           "count_plan","mean_premi27","d_mean_premi27","d_bench_premi27","premi27_bench","YEAR")

mfj.u <- rbind(mf15u,mf16u,mf17u)
mfj.nmkt <- rbind(mf15,mf16,mf17)
mfj <- rbind(mf15a,mf16a,mf17a)

### fit multilevel model
lout_d <- lmer(
  FAVOR ~as.factor(PID5)+as.factor(DATN)+EDUC+
    INCOME+BLACK+HISP+ASIAN+AGEN+COPCNHBL10+COPCHIS10+COMHINC10+COPCPOOR10+COPCUNEMP10+COPCUNEMPD1016+COPCBCH10+
    scale_log_count_plan+scale_mean_premi27+
    scale_d_mean_premi27 +d_bench_premi27+(1|YEAR),
  data=mfj %>% mutate(
    scale_log_count_plan=as.vector(scale(log(count_plan+1))),
    scale_mean_premi27=as.vector(scale(mean_premi27)),
    scale_d_mean_premi27=as.vector(scale(d_mean_premi27)),
    scale_d_bench_premi27=as.vector(scale(d_bench_premi27)),
    scale_premi27_bench=as.vector(scale(premi27_bench))
  )
)

stargazer(
  lout.nomkt,lout,lout.u,lout_d,
  digits=2,
  star.cutoffs = c(0.05, 0.01, 0.001),## ,star.cutoffs=,
  ## omit=paste(c(
  ##   "DATN","EDUC","INC","BLACK","HISP","ASIAN","AGEN","CO","PID",
  ##   "Constant"
  ## ), collapse="|"),
  ## covariate.labels=c(
  ##   "\\multirow{2}{4cm}{Number of plans (logged, in sd's)}",
  ##   "\\multirow{2}{4cm}{Mean Premium (in sd's)",
  ##   "\\multirow{2}{4cm}{Mean Change in Premium (in sd's)}"
  ## ),
  omit.stat=c("ser","f","adj.rsq", "rsq"),
  table.layout = "ts", float=F,
    out=paste0(
        the_prefix, "/figs/",
        "table3_tableA14_exchanges_multilevel_models.tex"
    )
)
